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A simple compact formulation is presented for the eigenvalue A of the Mathieu differential equa- 
tion y"-\-(k — 2h 2 cos 26)y = 0. The relations obtained thereby are employed to examine certain charac- 
teristic features of the solutions y — e il ' e V c„e 2i " e in both stable (v real) and unstable {v complex) regions. 

In the vicinity of the boundaries (v = m — integer) which separate these regions, the expressions derived 
for k are nonsingular and continuous; at the boundaries, these expressions are identical to available 
expansions which become therefore special cases of a more general formulation. The results are used 
to determine the magnitudes of the harmonic coefficients c„ and the instability intervals. The latter 
are shown to be intimately related to restrictions on the complex value of the characteristic exponent 
v. In addition, simple approximations that are suitable for perturbation problems (h 1 small) are derived 
for all of the characteristic quantities. 



1. Introduction 
Hill's differential equation and its special case, the Mathieu equation with real coefficients 

^ + U-2/* 2 cos2(9)y = ( L1 ) 

were first examined in connection with the vibration of elliptic membranes [1] ] and the motion of 
celestial bodies [2]. More recent studies were motivated by applications in the fields of acoustical 
and electromagnetic diffraction, frequency modulation, electromechanical transducer design, 
elliptic waveguides, etc.; since it is impractical to provide a complete bibliography within the 
present framework, the reader is referred to several excellent textbooks [3-6], which contain 
lengthy lists of references. 

Most of the investigations of the Mathieu equation (1.1) were concerned primarily with only 
a very limited range of solutions — the periodic ones. Thus, if the solution to (1.1) is expressed 
in the Floquet form 

y{e)=e ±[ ^ ]T c n (±iv)e 2ine (1.2) 

n~-oc 

where the ± signs relate to the two possible characteristic solutions, these will be periodic when- 
ever the characteristic exponent v is an integer (v = m = 0, 1, 2, 3, . . .). In those restricted 
cases, the solutions y{6) were extensively tabulated [7]. 

Of a somewhat different type are problems which do not require periodic solutions; in those 
cases, v is not necessarily a real integer and it may also become complex. Thus, the Mathieu 
equation (1.1) may express the propagation of electromagnetic waves in a medium that possesses 
dielectric or magnetic properties which vary periodically [8, 9, 10] ; in that case, 6 represents a 
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space variable and c n (Xiv) refer to the amplitudes of the space harmonics. An analogous situ- 
ation occurs in the case of certain types of oscillators [11] (particularly those involving parametric 
action) for which 9 is interpreted as a time variable whereas c n (— iv) relate to the amplitudes of 
the temporal harmonics. In all of these cases, tabulated results are few and restricted in range 
[6, 12]. 

The significance of the characteristic exponent v can usually be related directly to the wave 
character of the physical phenomena described by the Mathieu equation (1.1). Hence the deter- 
mination of v, i.e., its relation in terms of \ and h 2 op, alternatively, the determination of k=k(v, /i 2 ), 
are usually of considerable importance. It is also recognized that, when a forcing term is present, 
the source (generator) element usually prescribes the wavelike nature of the phenomena as re- 
flected through v\ hence \, rather than v, represents the unknown physical quantity to be deter- 
mined. One may then employ expressions that are found in the literature [1—7] and which relate 
k in terms of power series of h 2 for integral values of v. For v not an integer, the power series 

^ = ^ 2 + ^l)^ + 32(, 2 -iJi-4) /l8+ - • • d-3) 

is available; however, this relation converges poorly for v close to an integer and its applicability 
is therefore limited. 

The present paper presents a method for deriving explicit expressions for k = k(y, h 2 ) which 
hold for arbitrary values of v. The expansions involved are based on the continued fraction rela- 
tions [3-7] of the Mathieu equation and is an extension of a method for obtaining power series for 
v integer [13]. The expressions obtained herein are not singular for v integer, they are applicable 
for both real and complex values of v and appear to converge rapidly. A different method which 
also converges rapidly was recently proposed by Schafke, Ebert, and Groh [14]; in contrast to their 
results, the expressions presented here are explicit and they are particularly simple when approxi- 
mate expressions are required for small h 2 . 

In addition to expansions for the eigenvalue k, other characteristic parameters can be found 
via the expansions given here. Thus, the extent of the range for complex values of v is discussed, 
as well as the relation between all of the harmonic amplitudes c w (± iv). It is also shown that, in 
the limit of v becoming an integer, the results presented here are identical to those already avail- 
able [1—7] for the Mathieu functions and it is suggested that the present method leads to a more 
systematic approach for future computations of the pertinent series. 

Certain physical interpretations, notation and the representation of the various possible solu- 
tions are discussed in connection with the so-called "stability chart" in section 2. The derivation 
of the expansions for k = k(v, h 2 ) is presented in section 3 and the pertinent formulas for small h 2 
are given in section 4. Simplified expansions for v integer are discussed in section 5 and the range 
for v complex is evaluated in section 6. The magnitudes of the harmonics c n ( — iv) are dealt with 
in section 7. Finally, some comments are made in section 8 concerning the range for h 2 wherein 
the relations presented here are expected to be valid and to converge rapidly. 

2. Stability Chart 

The harmonic amplitudes c n = c n {±iv) are interrelated via the infinite continued fraction ex- 
pansions [3-6]: 

Cn h 2 \ h* \ h 4 \ 

Cn— 1 Lin \L*n+\ |-f>n+2 

c„ h 2 \ h 4 \ h 4 \ 
T~ = j-\j ~\j~ (2-2) 

C n +1 L/n \Ltn-l \Ltn-2 



where 



L n = (v+2n) 2 -k. (2.3) 
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By combining eqs (2.1) and (2.2), one obtains 



Ln = 1 1 
h 4 L n +i 



L n+2 



_LI. 



Ln- 



\L n - 



(2.4) 



which is essentially an equation for k = k(v, h 2 ). By solving it, one may obtain the eigenvalue X 
and thence the harmonic amplitudes c n via (2.1) or (2.2). 

Relation (2.4) is often presented [3-9] in the form of a so-called "stability chart", as shown in 
figure 1. This diagram illustrates some of the properties associated with the 4 eigenvalue X and the 
characteristic exponent v. Thus, for h 2 and X real, v may be either real or complex. In the latter 
case: 

v = m + ia (2.5) 

where a is real and m = 0, 1, 2, 3 . . . . For simplicity, the real part of v will always be taken as 
positive, in both the real and the complex v cases. 

Regions with real and complex values of v are called "stable" and "unstable" regions (or 
bands), respectively, in the literature [3-12]. A systematic definition of this terminology can be 
traced to van der Pol and Strutt [15]; the justification is afforded by the fact that, if v is complex,, 
the characteristic solutions in eq (1.2) are not bounded as 1 6>| — ^ °°- It is recognized that the 
"unstable" solutions should not necessarily be neglected in actual physical situations. Thus, when 
6 is bounded, both complex values of ±v are admissible; also, in certain unbounded situations, one 
of the solutions may be retained and only the other need be rejected. The latter case is illustrated 
by the radiation condition in electromagnetic fields [10]. 

From figure 1, one sees that unstable regions (bands) alternate with stable ones. The lines 
v = m consist of two separate branches (for all m^O) which join on the X axis and contain a region 
with v complex (=m-\-ia); it may be shown that the two v=m lines touch each other on the X-axis 
with an order of tangency of wi—l. As a consequence, for a fixed h 2 9 the width of the unstable 
band decreases as m increases. 

As indicated in figure 1, the curves for v = const, exhibit different variations in the stable and 
unstable regions. Those with v real start from the X-axis (at i> 2 = \) and proceed towards increas- 
ing values of h 2 \ on the other hand, the curves for v complex exhibit a minimum with respect to 
h 2 and are double-valued with respect to X. 

Whereas the stability diagram may be used directly [9, 10] to obtain one of the parameters 
X, v, or h 1 when the other two are given, this procedure usually leads to a graphical construction 
and therefore its accuracy and convenience are somewhat limited. Instead, simple analytical 
expressions for X = X(^, h 2 ) may be used and their derivation is given in the following section. 



Figure 1. The Mathieu stability diagram. 
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3. Relations for the Eigenvalue X 

In the present section, relations in the form \ = k{v, h 2 ) are obtained via the continued frac- 
tions of eq (2.4). When properly expanded as power series, these relations lead to equations for 
the eigenvalue \ which, unlike eq (1.3), are not singular at the band edges {v= m = 0, 1, 2, 3 . . .). 
In order to better appreciate the procedure employed here, it is pertinent to illustrate first the 
derivation of eq (1.3). 

By setting n = in eq (2.4), one obtains 

2 , *j M m .Jd _^J *_J 



which may be written as 



fhere 



k = v*-h\K + + K-) ,3.2) 



k , = _lJ_£L^„ 



1 Vj±2 K-»±3 



(3.3) 



Equation (3.2) is now viewed as an iteration formula 

K M = v 2 - h*[K+(kd + K-(kd] (3.4) 

for finding Xj+i when a trial value A;U' = 0, 1, 2, 3 . . . .) is assumed. To obtain an analytical result 
in terms of a power series in /i 2 , it is necessary to expand /£* as a power series. This is achieved 
via a useful series expansion for continued fractions [16] which, in the present case, is written as 

1 A* h» 

&±=— + — 






This expansion may be cast in the more compact form 

K± = ^B ±k B ±(k+l) h*« (3.6) 



k 



where B =l; B ±k =T[ A±j (3.7) 



j=i 



1 I /i 4 I h 4 \ h 4 \ h 4 

and ^^z^-^"^)"- ■ -~i^r\LZ (3.8) 

(/=1,2,3. . .;£ = (), 1, 2, . . .) 

and it is understood that A ± \=- — 

L,±\ 

It is recognized that the expansion in eqs (3.5) or (3.6) would be a power series in terms of h 4 

if terms with h 4 did not also appear in A±j. However, if h 4 is sufficiently small, one may employ 

the binomial series expansion to convert the above in a power series. Thus: 
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LxL\L* 



1 — 



+ 



L1L2/ L\L 2 L,i \ Lj\Li 






+ 1 + 



u) 



h* 



L\L\U 



+ 



L>L:s 



1- 



UL-. 



+ 



1+^+2-^-' 
Li L\L> 



L\L 4 

L\U) L*L*L*U 



0(/* 16 ) (3.9) 



and a similar power series for K~. Introducing these series into eq (3.4), one obtains the required 
iteration formula which yields X in terms of v and a power series in h 4 . The proper starting value 
for the iteration process is obviously X — v 2 for the stable regions and \o = m 2 (m = (), 1, 2 . . .) for 
the unstable regions, since these are the values of X at h = (see fig. 1). It is observed that by 
taking these starting values of X , the iteration formula (3.4) is in a form such that X, is obtained 
to 0(h 4i ). Hence, if X is to be found to 0[/i 4(i+1) ], it is necessary to find the series expansions for 
£± only up to 0(h 4i ). 

By using these considerations, one may verify that eq (3.2) yields formula (1.3). As already 
noted, however, this formula yields good results only for values of v which are not too close to a 
nonzero integer. When this condition is satisfied and if h 2 is not too large, both the iteration (3.4) 
and eq (1.3) converge rapidly. This fast convergence is due to the fact that X /+ i in eq (3.4) was 
isolated out from the term L , as indicated in eq (3.1) and, of all L n terms, L varies most rapidly 
with both h 2 and v, provided < real part of v < 1. This property is verified by noting that, for 
small /i 2 , one may write X = v 2 to yield 



L n = (v + 2nf-v 2 = <\n{v + n) 



(3.10) 



which is most sensitive for n = and n = —v = —m, where m is the integer that is closest to the real 
part of v. Hence, as long as v is real and sufficiently smaller than unity, L is the only sensitive 
term. However, in general both L and L- m may become critical; hence, an iteration formulation 
such as eq (3.4) would have to be based on a value for X;+i which is isolated out from both L and 
L- m . This is accomplished herein by employing certain properties of the pertinent functions L n , 
as discussed below. 

Consider first the problem of isolating the term Lq out of a continued fraction such as eqs 
(2.1) or (2.2). If \n\ is large, L will appear far down in the expansion and the process of finding 
its value will become quite tedious. It is therefore convenient to extend the concept of K + of 
eq (3.3) and define 









\L n +2 



(3.11) 



where the subscripts of the denominators form an increasing sequence. The term L n will appear 
only if n «£ 0; one may then write (see appendix): 



B 2 h~ 4n 
K n = A n + c " L4 (n«0) 



S-C„h* 
where the terms A„ and B„ are given in eqs (3.8) and (3.7), respectively, with 



(3.12) 



S = L -K t h* 

-(n+l) 

^o = C = 0; C n = YB- k B- (k+l) h 4k (n<0) 



(3.13) 
(3.14) 



and it is noted that the terms A n =A n (v), B n = B n (v) and C n = C n {v) were defined in such a manner 
that, if expanded as power series in h 4 , their expansion will start with a free term. 
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In addition to the simplification realized by writing K n in the form of eq (3.12), the terms 
L n possess a certain symmetry which is basic in the derivations obtained here. This feature is 
exhibited by taking v as 

v = m + A (3.15) 

where 

(3 (for v real) 
A= (3.16) 

ia (for v complex). 

For complex values of v, eqs (3.15) and (3.16) correspond to eq (2.5) and therefore m denotes 
the pertinent unstable region. In the stable regions, m is chosen so that |/3| ^ — ; hence m is then 
the integer closest to v in an arithmetic sense. 
Defining the "complementary" value of v as 

v = m — A = 2m — v, 

one recognizes that v corresponds to the complex conjugate value of v if the latter is complex. 
Consider now any function F = F(v); then 

F=F(v) = F(v), (3.18) 

where = indicates equality by definition. If F(v) is a polynomial, one may separate it as 

F = Fe + Fo (3.19) 

where Fe = Fe{m, A) and Fo = Fo{m, A) are the functions containing only even or odd powers of 
both m and A. For complex values of v, Fe, and Fo refer to the real and imaginary parts of F, 
respectively. One then also has 

F = Fe-Fo (3.20) 

h = FF={Fe) 2 -(Fof (3 .21) 

and it is recognized that, if v is complex, the modulus F is identical to the absolute value \F\. 
Applying definitions (3.17) and (3.18) to L n in eq (2.3), one gets 



L n =L_ (m+n) (3.22) 



which is the "symmetry" feature alluded above. Its importance lies in the fact that, whereas 
only the first fraction (with increasing indices in L n ) of eq (2.4) is expressed via K n of eqs (3.11) 
and (3.12), the second (decreasing) continued fraction is expressible in terms of K n via eq (3.22). 

A crucial point in the above definitions is the fact that X is regarded as an independent param- 
eter when taking the complementary of F(v). In particular, one notes that L n = L n (v, \) which may 
also be written_as L n = L n [v, \(v)]; definition (3.18), however implies that L n = L n (^, k) = L n [v, X(v)] 
and therefore L n # L n \v, X(v)]. 

The following discussion is more convenient if even and odd values of m are viewed separately. 
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(a)modd:m = 2r— l(r=l,2,3 . . .) 

For odd values of ra, eq (2.4) may be written, with n = 1 — r, as 



U- r 



U- r |UIr 



h A \ _ h x | h A | 

L r \T~Tr \TTr 



(3.23) 



which, by using relations (3.11), (3.18), and (3.22) becomes 



K x . r h\ 



£i-, 



(3.24) 



The sensitive L- m term is thus replaced with L . One recalls that, for an Iteration formula 
similar to that of eq (3.4), the term \,+i needs to be isolated from the L<> term, i.e., out of S of eq 
(3.13) in the present case; in particular, since now m 7^ 0, A;+i will be obtained from the even part of 
L . Writing 



X = S-Ci- r h 4 

the required iterated eigenvalue is given, via eqs (2.3), (3.13), and (3.25), by 



(3.25) 



(3.26) 



where k {m) will henceforth be used to denote the pertinent eigenvalue appropriate to the mth 
region, as defined in eqs (3.15) and (3.16). To find Xe in eq (3.26), it is realized that eq (3.24) leads, 
via eqs (3.12) and (3.25), to 



B 2 h 4(, '- l) 



X 



X 



(3-27) 



and the last relation is recognized to be a quadratic equation for Xe. By solving it and substituting 
into eq (3.26), one finally obtains: 



\< ± 2'-i>=(2r-l) 2 + A 2 -G'e-A 4 ±\[F' • A**'-"]** [2(2r- 1)A- Co • ft 4 ] 2 



(3.28) 



/here 



F' 



B* 



l-Al r h< 



G f = K { +C X -r-A^ r F'h«> 



(3.29) 



Result (3.28) is the required iteration formula for m — 2r— 1; though more complicated, this 
result is analogous to the iteration in eq (3.4) for m = 0. The ± signs in X. (, J ) refer to the two admis- 
sible eigenvalues for given h 2 and A. As already discussed in conjunction with the stability chart, 
all k ( '" ] (m/0) are double-valued for a given v and h 2 in the unstable regions; the + and — signs 
then refer to the larger and the smaller values of k {m \ respectively. In the stable bands, k {m) is 
single-valued and the two eigenvalues in eq (3.28) arise due to the choice of m in eq (3.15); it is 
then easy to verify that A.^"* correspond to /3 ^ 0, i.e., X ( ^ w) applies if m < v < ra + ^, whereas k (m) 

. . r 1 

is pertinent 11 m — - < v < m. 
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(b) m even: m = 2r (r — 1, 2, 3 . . . .) 

For even values of m, eq (2.4) is written, with n = — r and relations (3.11), (3.18), and (3.22), as 

L-r=(K l - r + K l - r )h* (3.30) 

which, by employing the notation of eqs (3.12) and (3.26), becomes 

IB 2 Z? 2 \ 
L- r -2Ae 1 - r h 4 =[^+^j h* (3.31) 

One then obtains again a quadratic with respect to Xe which yields a similar, but different, solu- 
tion when compared with the m odd case. By noting the value of L- r , as defined via eq (2.3), one 
now obtains 



\( ± 2r) = 4r 2 + A 2_ G '' e /i 4 ± J(/W A 4r)i + (4rA-G"o • h 4 ) 2 (332) 



where 

F" 



B 2 



tf-Xm-Uet-rh 4 (3.33) 

G , = Ki + Ci- r + F'h« r - 1 \ ( 3 >34) 

Equation (3.32) is the required iteration formula for even values of m; the significance of the 
± signs in X. { i n) was already discussed in the paragraph following eq (3.29). 

4. Formulas for the Eigenvalue X (w) for Small h 2 

Relations (3.28) and (3.32) are directly applicable for obtaining X (m) to any desired accuracy 
by means of numerical methods. In addition, analytical formulations may be obtained for X (m) 
provided h 2 is sufficiently small, this being particularly significant in certain physical situations. 
Thus, small h 2 occurs in problems wherein the periodicity with respect to the independent variable 
6 appears as a perturbation in an otherwise uniform configuration (A 2 = 0); it is then clear (e.g., 
from the stability chart) that instability phenomena may appear in certain ranges even if the 
perturbation h 2 is exceedingly small. 

To examine these effects, all of the functions appearing in eqs (3.28) and (3.32) are expanded 
as power series via relations which are analogous to eq (3.6) and (3.9). Using considerations simi- 
lar to those described in the paragraph following eq (3.9), formulas (3.28) and (3.32) are then em- 
ployed to yield relations for X (m) up to any desired power of h 2 . This method is applied in the pres- 
ent section to obtain expansion for A. (m) for small h 2 up to and including 0(h 4 ). It is easy to verify 
that, if all terms of 0(h 6 ) or smaller are neglected, it is sufficient to retain only first terms (h°) in 
the expansions for K u A n , etc.; the approximations involved are commented upon in section 8. 
The following results are thus obtained: 

(a) In the stable regions (v real). 

Formulas (3.4), (3.28) and (3.32) are used and expanded to 0{h 4 ); the iterations need be em- 
ployed only once and with a trial value of Xo (m) — ^ 2 - The exact definition of L n is important in these 
derivations, as commented upon in the two paragraphs following eq (3.22). 

- = 0: O^f; ^-- 2 ~2(I^ + 0(t^^)^ W-D 

1 1 ^ 3 

m=l: -=^_^1, l^p+^~; 



W=l + (l-v ± ) 2 



(5-v ± )h 4 



16(1 + i/ ± )(2-i/ ± ) 



2 + 



3h 4 



2 /I.6 



(1-^) 2 ±0(^I, (4.2) 



16(l + ^ ± )(2-i^ ± )J v ~ \64/' 
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3 n n 5 



X< ± 2) = 4+(2-^) 2 -h 



(7-^)/i 4 
(l-^)(3-^) 



- / h 8 
= V 16 ( 1 — ^ 



(l-^±) 



:+ 



4 + 



6(l + ^)(3-i/ ± ). 



^/379±384 L1 



(4.3) 



1 , ] 

m^ 3: m — —^v-^m,m^v + ^m-\- — 



X(«)=y2 



2(1-^1) 



+0 



5/n 2 + 7 



32(™ 2 -l) 3 (™ 2 -4) 



(4.4) 



It is noted that, to 0(/i 4 ), only the results for m— 1 and m = 2 are different from eq (1.3). On 
the other hand, eqs (4.2) and (4.3) imply that eq (1.3) may lead to very erroneous results even for 
small perturbations; this relates to the regions which are close to v=l and *> = 2, respectively. 

The order of magnitude of the first neglected term in the above expansions is also indicated; 
the figures give the exact result for v = m and it is expected that they yield a very good approxi- 
mation for v ¥" m. 

(b) In the unstable regions (y = m + ia). 

As illustrated in figure 1, h 2 and v may not be prescribed in a completely arbitrary fashion 
since, for m#0 and a given hr, a cannot exceed a certain value if X (m) is to be real. The limiting 
a corresponds to the minimum of the constant-^ curves in the unstable region; this value of a 
may be calculated by stipulating that the radicands in eqs (3.28) and (3.32) vanish, thus obtaining 



a *£ w m h 2m + 0{h 2m+2 ) 



(4.5) 



where the equal sign refers to the limiting value and 

1 



2m[2"'~ 1 (™-l)!] 2 



(for all m * 0). 



(4.6) 



Since the coefficients w m form a descending sequence, a is restricted to exceedingly small values 
for large m. This behavior is directly related to the extent of instability since a refers to the ex- 
ponential increase of the solution (1.2); consequently, large values of m (and therefore of \ (m) ) can 
accommodate only "weak" instabilities (at least for sufficiently small values of h 2 ) . 

With restriction (4.5) being accounted for, X (m) may be obtained in the same fashion as in the 
stable regions, except that the trial values for the iteration are taken as X ( [j ) = — a 2 {m = 0) and X ( J ) 
= m 2 (m 9^ 0). It is then obtained: 



-0: 



X (0) = -a 2 - 



2(l + a 2 ) 



+ 



(7- 5a 2 ) A 8 



32(1 -ha 2 ) 3 (4 -ha 2 ) 



(4.7) 



m= I: 



XL 1 



(S-a 2 )h 4 



(4-ha 2 )(16-ha 2 ) 



A 4 - (2a) 2 



1+- 



3h 4 



4-ha 2 )(16-ha 2 ) 



-O 



64 



(4.8) 



-2: 



x(2)-a » I 2(12 + a 2 )^ // A* \ 2 " 

^ ) - 4 -- 2 + ( 4-h a2)( 36 + a 2 ) ± Vl4Ta^j- (4a) 



1-h 



2h 4 



(4 -ha 2 ) (36 -h a 2 ) . 
379±384 IU 



-0 



13824 



m^3: Xl m > = m 2 -h 



2(m 2 -l 



r+O 



5m 2 -h 7 



32(™ 2 -l) 3 (m 2 -4) 
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A 8 



(4.9) 



(4.10) 



Due to analogous considerations, the comments following eq (4.4) (with respect to m ^ 3) also 
apply here. 

5. Characteristic Curves for the Periodic Mathieu Functions (v = m) 

As a particular case, formulas (3.28) and (3.32) are applicable for obtaining the characteristic 
curves of the periodic Mathieu functions {v = m). Observing that A = in these cases, all of the 
odd functions Ao n , Bo n , etc., are identically zero and the formulas simplify considerably. 

The notation for the band edges is generally given as a m = \ ( ™\v) \ v = m and b m = \. im) (v) | v = m in 
the literature [3-6] . One then obtains o 

ao = -2^ 1 /i 4 , (5.1) 

D2 /,2(2r-l) 

b- lr - l = {2r-\Y-{K,+C i -rW- °'-^ ^ , (5.2) 

#2 /,2(2r-l) 

a 2r _, = (2r- l) 2 - (Kt + Ci-r)h* + "^ ^ , (5.3) 

b ir = W-{Ki + Ci-r)h\ (5.4) 

a ir = W-(K 1 + C 1 - r )h*+ , D ]r^ {5 ^ 
(r=l,2, 3. . .). 

These relations may be used to generate the Mathieu function expansions for the eigenvalues 
a m or b m in terms of a power series in h 2 . It is then evident that certain known features are already 
apparent in formulas (5.1)-(5.5); thus, odd values of m lead to expansions in terms of powers of h 2 
whereas even values of m lead to expansions in terms of powers of h 4 . In fact, if these expansions 
are obtained iteratively with a first trial value of X = ^ 2 , the above formulas yield results which 
are identical with expressions available in the literature [3-7]. The identification follows upon 
recognizing that, for example, eqs (5.2)-(5.5) are a restatement of eqs (20—2, 21-24), respectively, 
on p. 723 of reference 6. 

The advantage of the present formulation is evident in that it represents a systematic technique 
for generating the power expansions for k {m) . Whereas it is recognized that all methods (including 
the present one) become quite involved and tedious if a m or b m need to be found in terms of large 
powers of A 2 , it is nevertheless clear that formulas (5.1)-(5.5) are already expressed in such a 
concise and systematic manner that considerable simplification occurs. 

6. Instability Intervals 

The width of the band between each pair of the characteristic v = m curves is of special im- 
portance since it yields information on the extent of the range of instability. For an arbitrary 
value of h 2 , this instability interval is defined by 

and is easily determined via relations (5.1) to (5.5) which yield: 

2B 2 

I— A 2 h 4 ^ 6 *^ 

1 - r 

no 



Both results may be expressed also as 



W m = 2 



|_2»<- 1 (m-l)! 



0(/* 2 ™+ 2 ) (6.4) 



which is obtained by taking the first term in a series expansion of eqs (6.2) or (6.3) in powers of h 2 . 
It is observed that, to this order of approximation 

W m = 4>mw m h 2m (6.5) 

where w m are the coefficients in eq (4.6). Result (6.4) is identical to that already obtained by Bell 
f 17], Levy and Keller [18] and Hochstadt [19] who used different approaches. In the present case, 
the instability interval W m turns out as a simple corollary to formulas (5.2)-(5.5) and the accuracy 
of eq (6.4) can be estimated as discussed in section 8. It is also noted that W m is simply related to 
the maximum allowable value of a since, introducing eq (4.5) into eq (6.5), one has 

W 

«^if + <W 2M + 2 )- (6.6) 

Hence, as m increase, the range of a becomes more restricted than that of W m . Phrased dif- 
ferently, larger values of v (or k) lead to regions with decreasing unstable intervals; the extent of 
instability itself is proportional to the interval W m , but inversely proportional to v (or k). 

7. Magnitude of Harmonics 

The magnitude of the harmonic coefficients cj±iv) of eq (1.2) may be found by using the 
appropriate values of \ (m) obtained above. For simplicity, the + sign in c n (±iu) will be implied 
subsequently; hence, c„ is taken to denote c n (iv). To find the coefficients c»( — iv), one may employ 
the identity c n ( — iv) = C- n (iv) which is proven in the literature [3—6]. 

Using eqs (2.1), (2.2), (3.11), and (3.22), one gets 

7~ = ~ h2Kn (7.1) 

Cn C (m + n) » 2 ~77 

C n +l C-(m+n+l) w '^' 

As clarified below, the c ( > and c_ m harmonics turn out to be particularly important. For these har- 
monics, one obtains the following relations: 

(a) For m odd; m = 2r-l(r=l, 2, 3 . . .). 
By using eq (7.1) successively, one obtains 

-H rr7=(-* , )- , n^ (7.3) 



and eq (7.2) similarly yields: 



so that, by division: 



C] 



-=n7r= { - h2y n K > (7.4) 



o, = C\-2r_ 



(-h 2 K (ior r=l) 



r-2 K_. 



757-615 0-65-8 



-i^nf(f<>rr>l) ,7.5) 
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(b) For m even; m = 2r(r= 1, 2, 3 . . .). 

In a manner similar to that for m odd, one gets: 



C-2r . 
C-r 



2r ~ . r-1 _ 



(7.6) 

(7.7) 



1=0 



and, by division: 



*-— -Ifc- 



Co 



(7.8) 



By using eqs (7.1), (7.2), (7.5), and (7.8), one may now express all c n by means of the following 
relations: 



(~h*)*l[Ki (n>0) 

m- 1 

(-/i 2 )-" ] [ K-i (-r<ra<0) 

i=m+n 

m—l 

(- /i 2 )™ +w # J [ K-i (-m<n^-r) 

i=—n 

-(m+n) _ 
(_ A 2)-(m+n)/ ? ^. ( n <--|») 



Co 



(7.9) 



and it is understood that R stands for R' or R" for odd or even values of m, respectively. The 
fundamental harmonic c may be taken as unity or normalized via an orthogonal relationship 
[3 — 6]; all these normalization procedures, however, lead to a real value for Co- 
Equations (7.5) and (7.8) lead to the useful relation 



R2 = RR = l 



(7.10) 



which is obvious for R = R"; for /? = /?', this is verified by observing that eq (3.24) accounts for 
the extra K\- r term in eq (7.5). One therefore gets 



r 2 = ~r 2 / 

c n c -(m+n) 



(7.11) 



which is obtained via eq (7.9) and (7.10). 

Result (7.11) is interesting since one recalls that, in the unstable bands, any function F(v) 
leads to a complementary function F{y) = F{v) which is identical to F*{v) =F{v*), where the 
asterisk denotes the complex conjugate. Consequently, the modulus F(v) is then identical to the 
absolute value \F(v) |, so that 



\C n I — |C-(m+n)| 



(7.12) 



in the unstable bands. The harmonic amplitudes are then equal in pairs, as shown in figure 2. 
Relation (7.12) holds up to and including the band edges {v = m) where all variables are real; hence 



C n — — C-(m+n) 
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(7.13) 



since all terms are real and therefore any function F(v) =F*(v) =F{v) at the band edges. The 
proper sign in eq (7.13) is determined by the sign of R =± 1 itself at the band edges. For m odd, 
the sign of R' is given by the sign of — h 2 Ki- r = — h 2 K\- r because the product term in eq (7.5) is 
obviously + 1 at the band edge. By using eqs (3.12), (3.13), and (3.24), one has 



-h 2 K^r=-A^- r h 2 



B 2 /i 2(2r " 



(2r-l) 2 -^ 2r - 1) -(K l -i-C 1 -r)h 4 



R' 



1 



(7.14) 



where X^ 2r ~^ refers to bzr-i and X ( + r_1) refers to «2r-i of eqs (5.2) and (5.3), respectively. It is easily 
seen by inspection that, in order to comply with eqs (5.2) and (5.3), R' =+ 1 need be taken for ct2r-\ 
whereas /?'=— 1 is appropriate for b 2r -i in eq (7.14). 

The even m case is somewhat more complex. Checking for K\- r , it is seen that ai r of eq (5.5) 



leads to a finite value for K\- r but bzr of eq (5.4) requires that \K\ 



►oo. Hence R" = + 1 for a> r 



(since then all K-t=K-i are finite) but is yet undetermined for 6 2 r- To clarify the latter situation, 
one uses eq (3.30) to find the limit 



im 



Kt- 



A, 



A^O^l-r 



Tim 



L- r 

h A K X - 



1 



1 



and since all K-i = K-i (i ^ r— 1) in eq (7.8) at the edge pertinent to 6 2r , eq (7.15) leads to R" = — 1 
at that edge. In conclusion, one obtains that 



c n = J rc-( m+ n) for a n 
C w =— C-(m+n) for b, h 



(7.15) 



and therefore the function y{6) of eq (1.2) is in the form of a Fourier cosine or sine expansion, 
respectively, on the right or left edge of the unstable regions. 

A special feature of the b lr edge is that, since /Ci_ r — **>, one gets from eq (7.6) that c- r = which 
also agrees with the second of eqs (7.15); hence, the c~ r harmonic vanishes at the left edge of all 
the unstable regions which are characterized by even values of m = 2r. It is appreciated that 
this harmonic is exceptional in that it is the only one which (for m even) does not pair off with any 
other harmonic in the unstable regions but, as shown in figure 2 and eq (7.12), it can be said to 
pair with itself. 

-2rH o 



-2r+2 
■2r+3 



-r-2 -r+2 



II I l___ 

(a)m ODD (m = 2r-l) 



Figure 2. The magnitudes of the harmonic coefficients 
c n (iv) in the unstable regions. 
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To conclude the discussion, it is pertinent to find the magnitudes of the harmonics for small 
h 2 . In that case, one may approximate K\ in the unstable regions as 



^r = Li= (m + A + 2i) 2 - (m + A) 2 = 4i(m + i) = = - 



(7.16) 



For K-i, however, condition (4.5) precludes an easy approximate evaluation. The order of magni- 
tude may nevertheless be estimated by replacing i with — i in eq (7.16). One then finds: 



h 2 V 



4/ nl(m + n)l 
( 



0(h 2 ^ 



r° 



h 2 







m + n-l)\ 1 
-it)!(i»-l)!J 

(-n-l)l 



4 

^2\-(m+n) 



(m + ;i)!(m-l)!. 
' (-m-ji)I(-ii)!" 



(/i>0) 

(-r<n<0) 

( — m < n ^ — r) 
■ O (h 2 - m -») (n * 



(7.17) 



• m) 



for the unstable regions. Numerical and approximate analytical considerations indicate that, for 
m > 3, the order of approximation for — m < n < in eqs (7.17) approaches the actual value very 
closely. 

In the stable regions, property (7.10) does not hold any longer and K- ( turns out to be quantita- 
tively dissimilar from K-i. To find the latter, one notes that eqs (3.1) and (3.13) imply that 



S - C-ih 4 = U -K,h 4 - C-ih* = {K- - C-dh 4 



(7.18) 



where K~ and K + = Ki were defined in eq (3.3). Employing the definition of C-% in eq (3.14) and 
the series expansion (3.6) for K~, one finds 



K--C-i = B-iB- ii + 1) h 4i + 0[h« i +»]. 
Introducing eqs (7.18) and (7.19) into eq (3.12), one obtains via eqs (3.7) and (3.8) 
K-i=A-i + -z ri+0(A ) = ~ -rz \-0(h°) = ri \-0(h 2 ) 



(7.19) 



(7.20) 



In view of the comments following eq (3.22), K-\ cannot be obtained directly from eq (7.20) and one 
needs to use definition (3.12) instead. It is then observed that 



S = Lo-K l h 4 = (2m-p) 2 -\-K l h 4 ; 



(7.21) 



also, an inspection of eqs (4.1)— (4.4) for A reveals that S is of 0(A°), provided v is not too close to 
m (band-edge). Thus 



K-i=A-i + 0(h 4i )==^-+0{h 4 ) = — : 

L-j 4(ra — i)(m — v — i) 



+ 0(h 2 ) (i ^ m) 



(7.22) 



and this approximation is valid for use in eqs (7.9) since i ¥^ m therein. 

The terms Ki and K { are obtained in a similar form to that of the unstable regions. Thus, 



K t =Y +0(h*)= } , + 0(/t*), 

Lj 4t(t + v) 



(7.23) 
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^S +w " 4o»+o(L- y+ a +o»')- 



(7.24) 



Introducing the relations above into eqs (7.9) one has for the stable regions: 



c n ( h 2 \\»\ r(l± y) 



(7.25) 



where T(^) is the gamma function of argument x. 

When comparing the last results with eqs (7.17), one recognizes that the order in h 2 for |c„/c | 
is the same in both stable and unstable regions for n &* — r; for n < — r, on the other hand, |c„/c | 
changes from 0(h 2 \ m+n \) to 0(h 2 ^), i.e., by a factor of h 2m . This is due to the fact that pairing of 
the harmonics occurs in the unstable region only. 
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FIGURE 3. Variation of the harmonic coefficients c n {iu) 
within the stable regions, for h 2 =l. 
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The above considerations indicate the following: (a) for n^ — r, \c n \ decreases monotonically 
relative to c as k increases; (b) every harmonic \c n \ with n< — r increases near and at the band 
edges and decreases towards the middle of the stable bands. These remarks apply to cases with 
h 2 small and, as one follows the value of c n with n < 0, it is clear that this harmonic will behave in 
the manner (b) up to a value of X = r and in the manner (a) for X > r. For n > 0, the harmonic am- 
plitudes behave only in the manner (a) since then n is always greater than — r. 

The above features are illustrated in figure 3 wherein a few harmonic amplitudes are plotted 
for the range ^ v ^ 6 within the stable regions. The vanishing of the c~ r harmonic at the left 
band edge with m even is shown in the figure by c~i, c-2, and Cs which vanish respectively, at 
v = 2, 4 and 6. Also, relations (7.16) are verified; thus, |c-i|=co at v=\, |e_ ^ | = c at v = 2, \c\\ 
= |c_ 2 1 at v = 3, etc. The constant monotonic behavior of c n for n>0 should also be observed. 

8. A Note on Accuracy and Convergence 

In the foregoing derivations, the exact meaning of "/i 2 small" was not explored and this ques- 
tion is examined in the present section. 

It is first remarked that, if a numerical computation is desired, the results (3.28) and (3.32) 
may be valid for the calculation of X for any values of v and h 2 . One notes that all the expressions 
involved are in the form of finite continued fractions, except for Ki=K + which is infinite. The 
latter, however, is expandable via eq (3.6) into a series consisting of terms Kk that, for large k, 
vary as 



Kk = B k B k+1 h*x -* — — k » - 



4(£+l)(m + k+ 1) 
L k+1 || Lj 



kl(m + k)\. 



'(I 



(8.1) 



which converges very rapidly. Hence, for numerical computation, k is obtainable within any 
desired accuracy provided the iteration process itself converges. However, the convergence of 
the iteration process is still an open question. 

With the above observations in mind, the accuracy of the results for small h 2 can be examined. 
As already noted, only the leading terms in the expansions for K u Ai- r , etc., need be retained for 
results up to and including a term of 0(/i 4 ). With respect to K u this means that one needs 

h 4 <^ \LiLi+i | for all i > 0, (8.2) 

in order to justify the assumed approximation. All of the other functions Ai- r , fii-r, and C\- r are 
essentially expressible in terms of A\- T \ retaining only the first term in the latter means that 

A 4 <^ \L- } L X -j\ for all Kj < r (8 . 3 ) 

If one replaces X by the first approximation \ — v l = m 2 in the two inequalities, one gets: 

h 4 <^ I6i(i + 1) (m + i) (m + i + 1) {i > 0) (8.4) 

and h 4 < I6j(j- 1) (m-j) (m-j+ 1) (Kj < r). (8.5) 

The worst case for inequality (8.4) is m = 0, i = 0, i.e., 

h 4 < 64 (8.6) 

whereas, for (8.5), the worst restriction occurs at m = S and7 = 2 since only 1 <j < r applies; hence, 

h 4 < 12 X 32. (8.7) 
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FIGURE 1. Percentage error incurred in the value of A, for \i 2 — I. 

Solid lines refer to the error produced by eqs (4.1H4.4) and dashed lines refer to the error due to the 
(1.3). The dot-dashed curve relates to the error which occurs when eq (4.2) is used instead ofeq (4:1). 



ofeq 



Equations (8.6) and (8.7) indicate that the results obtained should hold well up to and including 
values of h 2 of order unity. An analogous examination of L (instead of L) leads to the same restric- 
tion. To check these considerations, the exact results for h 2 = 1 are compared in figure 4 with 
those obtained via eqs (4.1)-(4.4). The solid lines indicate the error produced by using eqs (4.1)- 
(4.4) whereas the dashed lines refer to the results obtained via eq (1.3). It is noted that consid- 
erable improvement is obtained, especially near v=l and v=2. It is also interesting to observe 
that smaller errors are obtained if eq (4.2) is used to calculate \ for ^ v ^ 0.5 instead of eq (4.1); 
this feature is shown in figure 4 by the dot-dashed line. Calculations were also carried out for 
h 2 = 0.1 and 2.5; the errors incurred then are respectively, less than 0.1 percent and 20 percent for 
all v, provided eqs (4.2)-(4.4) are employed (i.e., eq (4.2) is extended to cover the range < v < 3/2). 
The above comments indicate that the formulation discussed in the present paper is useful 
from both analytical and computational aspects. While it is recognized that the iteration process 
may not always be convergent, it is also probable that the present formulation might provide a 
tool for examining the range of validity for some of the expressions of the Mathieu functions and 
solutions whose convergence is still an open question. 
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9. Appendix 

Assume that K n , as defined in eq (3.11), may be written as 

Qn 



K n = P n + 



S-R„ 
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(n«0) 



(Al) 



where P n , Q n , and R n are functions to be determined in such a manner that the above becomes an 
identity; S was defined in eq (3.13). One then has 



K n -\ '■ 



1 



Ln-i — h 4 K n 



= Pn-i + 



Qn-l 
S — R n -1 



Substituting for K n , eq (A2) yields after some algebraic manipulation, 

On 



Pn-1 + 



1 



Qn-l „ 
S — Rn-1 L„-i — h 4 P n 



1+- 



Ln-l-h*P n 



S-R„- 



Qn 



Ln-l-h*P n 



(A2) 



(A3) 



and, identifying terms, one obtains 



p»-i= 

<?n-l = 





1 


L/n— 


,-MP» 




Qnk 4 


(L„. 


-l-h*P n f 


p 


Qnh* 



Ln-h*Pn 



(A4) 
(A5) 
(A6) 



The last relations form a system of difference equations which are solved by noting that the bound- 
ary condition at n = stipulates that 



Po=i?o = Oand<?o=l. 



Applying this to (A4) — (A6), one obtains easily 

P = JJ_^U_^1 
n t if \r^ 



h 4 \ h 4 



i=l 
R n = h 4 ^P.jQ:.j 



(A7) 

(A8) 
(A9) 

(A10) 



all of which hold for n < and it is understood that P-i= 1/L-i. Upon identifying 

An - Pn, B* - Qnh 4n and C n = Rnh-\ (All) 

it is readily verified that eq (3.12) is proven and that eqs (3.7), (3.8), and (3.14) are satisfied. 
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